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ABSTRACT 



A finite element formulation with linear triangular elements was used to 
solve the steady-state, two-dimensional conduction heat transfer equation 
in the condenser wall section of an internally finned rotating heat pipe. 

A FORTRAN program using this method was coupled with the ADS program 
for automated design of the internal heat pipe fin geometry to optimize 
heat transfer. An increase in surface area, which increases heat transfer, 
also increases the condensate level, which decreases heat transfer. 

The additional condensate level does not offset the advantage gained by 
the increased surface area. The investigation provided combinations of 
fin half angle, number of fins, and fin height for an optimum design. 

Water is used as the working fluid and the heat pipe is constructed from 



copper. 
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I. INTRODUCTION 



A. THE ROTATING HEAT PIPE 

The rotating heat pipe is a closed container designed to transfer a 
large amount of heat in rotating machinery. Since the heat pipe operates 
on a closed two-phase cycle, the heat transfer capacity is greater than for 
solid conductors. Also, the thermal response time is less than with solid 
conductors. The three major elemental parts of the rotating heat pipe are: 
a cylindrical evaporator, a truncated cone condenser, and a fixed amount 
of working fluid as shown in figure 1. 

An annulus is formed by the working fluid in the evaporator. This 
occurs at rotationary speeds above the critical speed. The addition of heat 
to the* evaporator vaporizes the working fluid. A pressure differential 
between the evaporator and the condenser causes the vapor to flow 
towards the condenser. The vapor is transported to the condenser with its 
latent heat of vaporization. Condensation of the vapor on the inner wall is 
caused by external cooling. This condensation releases the latent heat of 
evaporation. This condensate is forced to flow back to the evaporator by 
a component, acting along the condenser wall, of the centrifugal force 
which is caused by the rotation of the heat pipe. As the condensate 
collects in the evaporator the cycle is repeated. 

Since the evaporator and condenser portions of a heat pipe function 
independently, needing only common liquid and vapor streams, the area 
over which heat is introduced can differ in size and shape from the area 
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gure 1. Schematic Drawing of a Rotating Heat Pipe 
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over which it is rejected, provided that the rate at which the liquid is 
vaporized does not exceed the rate at which it can be condensed. 
Therefore, high heat fluxes generated over relatively small areas can be 
dissipated over larger areas with reduced heat fluxes; allowing a 
cylindrical evaporator and a truncated cone condenser. 

Capillary action acts to drive the condensate back to the evaporator 
in a conventional heat pipe. No limitation due to capillary action is 
encountered in a rotating heat pipe nor are external pumps or gravity 
depended on for the flow of the working fluid. Therefore, the rotating heat 
pipe can be used in any orientation [Ref. 1]. 

B. OPERATING LIMITS OF A ROTATING HEAT PIPE 

The first theoretical investigation of the rotating heat pipe conducted 
at the Naval Postgraduate School was performed by Ballback [Ref. 2] in 
1969. Various fluid dynamic mechanisms limit the performance of a rotating 
heat pipe. Ballback [Ref. 2] studied these mechanisms and an estimation of 
the sonic limit, boiling limit, entrainment limit, and condensing limit of 
performance was made. 

1. The Sonic Limit 

The maximum flow of the vapor is set by the choked flow 
condition in the rotating heat pipe. This limiting vapor flow rate occurs 
when the heat flux is increased and limits the amount of energy the vapor 
can transport. The rotating heat pipe effectiveness is reduced due to this 
limitation. The limiting heat transfer rate due to this condition is: 
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( 1 ) 



Q, ' P /M** 



and the vapor velocity is considered to be sonic. 



where 

U v = velocity of the vapor in ft/sec, and 

A = cross sectional area for the vapor flow in ft 

c = sonic velocity in ft/ sec 

g 0 = gravitational constant 

k = ratio of specific heats 

R = gas constant in ft-lbf/lbm R 

T = absolute temperature in degrees Rankine 

p v = density of vapor in lbm/ft 

2. The Boiling Limit 



The transition from nucleate to film boiling was hypothesized by 



Kutateladze [Ref. 3] to be a completely hydrodynamic process. He 
determined the following theoretical formula for predicting the burnout 
flux: 



1 



U v -c - (g 0 kRT ) 2 



( 2 ) 



1 




(3) 



where 



K = constant value 



0 

Au = heat transfer area in the boiler in ft 
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hfg = latent heat vaporization in Btu/lbm 
a = surface tension in lbf/ft 

2 

g = acceleration of gravity in ft/hr 
pf = density of fluid in lbm/ft 3 
p v = density of vapor in lbm/ft 3 . 

A constant value for K in the range of 0.13 to 0.19 was suggested by the 
experimental data obtained by Kutateladze. 

3. The Entrainment Limit 

The flooding constraint in a wickless heat pipe was examined by 
Sakhuja [Ref. 4]. The correlation he developed is: 

A x C 2 hJgD( P/ - Pv ) Pv 

V, J (4) 

L+(p J P )*f 

where 

Qt = heat transfer rate in Btu/hr 

2 

A x = flow rate in ft 

C = dimensionless constant, 0.725 for tube with sharp edged flange 
hfg = latent heat of vaporization in Btu/lbm 
g = acceleration due to gravity in ft/hr 
D = inside diameter of heat pipe in ft 
pf = density of the fluid in lbm/ft 3 
p v = density of the vapor in lbm/ft 3 . 
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4. 



The Condensing Limit 

The condenser section of a rotating heat pipe was modeled as a 



truncated cone by Ballback [Ref. 2]. Using this model, the condensation 
limitation for a rotating heat pipe was determined by Ballback [Ref. 2]. He 
developed the following condensation limit: 



where 

Qj. = total heat transfer rate in Btu/hr 

kf = thermal conductivity of the condensate film in Btu/hr-ft-F 



0 ) = angular velocity in 1/hr 

hfg = latent heat of vaporization in Btu/lbm 

T s = saturation temperature in F 

T w = inside wall temperature in F 

Uf = viscosity of fluid in lbm/ft-hr 

<(> = half cone angle in degrees 

R 0 = minimum wall radius in ft 

L = length along the wall of the condenser in feet. 

The geometry and speed of the rotating heat pipe, and the physical 
properties of the working fluid comprise the condensing limit equation. 



shown in Table I, the amount of heat that can be transferred from the 
rotating heat pipe is limited by the condensing limit. However, the 




3 

pf = density of fluid in lbm/ft 



For a rotating heat pipe with the physical characteristics as 
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limitations imposed by the sonic limit, boiling limit, and entrainment limit 
may become important as the heat pipe geometry and operating conditions 
vary. 

TABLE I. ROTATING HEAT PIPE SPECIFICATIONS 



Length 


, r^mi- -■ — ■■ » — - 

9.000 inches 


Minimum Diameter 


1.55 inches 


Wall Thickness 


0.03125 inches 


Internal Half Angle 


1.000 degree 


| Rotating Speed 


3600 RPM j 



To enhance the heat transfer capacity of the rotating heat pipe. 



internally finned condensers have been used to raise the condensing limit 
line. Thinner films occur near the ridges of the fins while thicker films 
occur in the troughs. The thinner film on the ridges provides a lower 
thermal resistance to heat flow, while the thicker film in the trough 
provides a higher resistance. A compromise between the improvement on 
the ridges and the degradation in the troughs is necessary for an overall 
heat transfer improvement [Ref. 5]. 

C. ANALYSIS OF THE INTERNALLY FINNED ROTATING HEAT PIPE 

Schafer [Ref. 6] developed an analytical model for a heat pipe with 
a triangular fin profile (figure 2). This model was developed in order to 
raise the condensing limit by the addition of internal fins. An assumption 
of one-dimensional heat conduction through the wall and fin was made for 
Schafer's model. 

A two-dimensional heat conduction model using a Finite Element Method 
was developed for this same case by Corley [Ref. 7], A parabolic 
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temperature distribution along the fin surface was assumed by Corley [Ref. 
7]. A significant improvement in the predicted heat transfer performance 
was indicated by his results. By using the two-dimensional model an 
increase of approximately 75 percent in the heat transfer performance was 
seen over the results from the use of the one dimensional model. However, 
Corley [Ref. 7] noted that at the fin apex an error as great as 50 percent 
was possible which could result in the total heat transfer being in error 
as much as 15 percent. 

A modification was made to Corley's computer program by Tantrakul 
[Ref. 8]. In order to minimize the heat transfer error at the apex of the 
fin Tantrakul increased the number of finite elements used. His results 
with this modification converged with the results of Corley. 

Purnomo developed a linear triangular finite element model (figure 3) 
used in a two-dimensional Finite Element Method solution. Purnomo's [Ref. 
1] Finite Element Method program also worked and converged. To maximize 
the heat transfer from the rotating heat pipe the condenser geometry was 
varied. Using Purnomo's code parametric studies were conducted. However, 
the best geometry was not indicated in these studies. Purnomo's code was 
written to perform one analysis at a time. Davis [Ref. 9] modified 
Purnomo's code to allow for numerous analysis to be made using the 
optimization code COPES/CONMIN. Davis' Finite Element . Method code 
incorporating the optimization worked and converged, resulting in an 
optimum design for an internally finned rotating heat pipe. 
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I 




Figure 2. Internally Finned Condenser Geometry, Showing Fins, 
Troughs and Lines of Symmetry. 
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b • haight of tha fin 

• condansar half angle 
a » half of the fin width 
c® half of the trough width 

Figure 3. Condenser Geometry Considered with 40 Linear Triangular 
Finite Elements. 

D. THESIS OBJECTIVES 

The objectives of this thesis were: 



1. To modify Davis' [Ref. 9] computer program so that it is compatible 
with the ADS (Automated Design Synthesis) program [Ref. 10] and 

can be used for analysis and automated design of ro tatin g heat 
pipes. 
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2. To use the resulting program to obtain an optimum design for an 
internally finned rotating heat pipe to obtain experimental data to 
compare with the analytical results. 

3. To use the resulting program to obtain numerical results in place of 
data obtained from costly experimental operations. 



II. NUMERICAL OPTIMIZATION 



A. BACKGROUND 

The parameter that is minimized or maximized during the design 
process is called the design objective. The design objective is minimized 
or maximized by changing the design variables within the design constraint 
limitations. This process is called numerical optimization. An assortment of 
physical, aesthetic, economic and, on occasion, political limitations must be 
met by the design constraints for the design to be acceptable. For the 
optimization process to work, the design criteria must be described in 
numerical terms. This is not always easy. 

A computer program can be written to perform tedious and repetitive 
calculations necessary to optimize the problem once it is stated in 
numerical terms. For this reason, computer analysis is commonplace in most 
engineering organizations. For example, in heat transfer design the 
configuration, materials, and method of heat removal may be defined and 
a finite element analysis computer code is used to calculate temperatures, 
heat transfer rates, and other response quantities of interest. If any of 
these parameters are not within prescribed bounds, the engineer may 
change the method of cooling or other defined quantity and rerun the 
program. The engineer makes the actual design decisions, the computer 
code only provides the analysis of a proposed design. This is the commonly 
used approach which is called computer-aided design. 
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Analysis codes are commonly used for tradeoff studies. For example, 
an analysis code might be run on the distance a truck can go on a tank 
of fuel. For different loads, different distances are calculated which can 
be used in a range-payload study. 

Fully automated design is the logical next step to computer-aided 
design. The computer makes the actual design decisions or trade-off 
studies based on input criteria in fully automated design. Minimal 
information is requested from the operator during the actual design 
process. Numerical optimization offers numerous improvements over the 
traditional approach to design. These improvements include: time reduction 
in design decision making; a rational, directed design procedure; and the 
procedure is unbiased by intuition or experience. The probability of 
obtaining a non -traditional solution is thereby improved. Engineering 
intuition and experience are still necessary to decide if the design obtained 
is an improvement and feasible. 

B. AUTOMATED DESIGN SYNTHESIS (ADS) 

Vanderplaats [Ref. 10] developed a general purpose numerical 
optimization program containing a variety of algorithms, ADS. ADS is a 
FORTRAN program that optimizes a numerically defined objective function 
subject to a set of constraint limits. The solution of the problem is 
separated into three levels: 

1. Strategy - Optimization strategy such as Augmented Lagrange 
Multiplier method or Sequential Linear Programming. 

2. Optimizer - Actual algorithm to perform the optimization. 

3. One-Dimensional Search - Line search routine used by optimizer. 
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Flexibility to solve a wide variety of engineering design problems is 
given by the combinations of nine strategies, five optimizers, and eight 
one-dimensional search options. The following definitions are necessary to 
discuss the use of ADS: 



1. Design Variables - Those parameters which the optimization program 
is permitted to change within allowed bounds in order to improve 
the design. Design variables appear only on the right hand side of 
an equation and are continuous. 

2. Design Constraints - An inequality constraint requires that some 
function of the design variable(s) remain less than a specified value. 
Design constraints may be linear or nonlinear, implicit or explicit, 
but they must be continuous functions of the design variable. 

3. Objective Function - The parameter which is going to be minimized 
or maximized during the optimization process. The objective function 
may be linear or nonlinear, implicit or explicit, and must be a 
continuous function of the design variables. The objective function 
usually appear on the left side of an equation. 



C. PROGRAMMING GUIDELINES 

Any computer code developed for engineering analysis should be 
written in such a way that it is easily coupled to a general purpose 
optimization program such as ADS. Therefore, a general programming 
practice is outlined here which in no way inhibits the use of the computer 
program in its traditional role as an analytical tool, but allows for simple 
adaption to ADS. 

ADS is called by a user-supplied calling program. ADS does not call 
any user-supplied subroutines. Instead, ADS returns control to the calling 
program when function or gradient information is needed. The required 
information is evaluated and ADS is called again. This provides considerable 
flexibility in program organization and restart capabilities. Various internal 
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parameters are defined on the first call to ADS which work well for the 
"average" optimization task. However, it is often desirable to change these 
in order to gain maximum utility of the program. Figure 4 is the program 
flow diagram for the case where the user wishes to over-ride one or more 
internal parameters, such as scaling, convergence criteria, or maximum 
number of iterations. 

After initialization of basic parameters and arrays, the information 
parameter, INFO, is set to -2. ADS is then called to initialize all internal 
parameters and allocate storage space for internal arrays. Control is then 
returned to the user, at which point these parameters, for example 
convergence criteria, earn be overridden if desired. At this point, INFO will 
have a value of 1 and the user must evaluate the objective function, OBJ, 
and constraint functions. ADS is called again and the optimization proceeds. 
Since, in this case, the gradient calculation control, IGRAD, has a value of 
zero, all gradient information is calculated by finite difference within ADS. 
When INFO has a value of zero, optimization is complete. 
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BEGIN 



DIMENSION ARRAYS 
DEFINE BASIC VARIABLES 
IGRAD=0 (USE FINITE DIFFERENCE GRADIENTS) 

INFO = -2 
CALL ADS (INFO...) 

IF INFO = 0, EXIT. ERROR WAS DETECTED 
ELSE 

OVER-RIDE DEFAULT PARAMETERS IN ARRAYS WK AND IWK IF 

DESIRED 

CALL ADS (INFO...) 



NO 



YES 



INFO = 0 



EVALUATE OBJECTIVE EXIT OPTIMIZATION 

AND CONSTRAINTS IS COMPLETE 



Figure 4. Program Flow Logic: Over-Ride Default Parameters, Finite 
Difference Gradients [Ref. 10]. 
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m. FINITE ELEMENT SOLUTION 

A. REVIEW OF THE PREVIOUS ANALYSIS 

As stated previously, the heat transfer solution for a one-dimensional 
model of an internally finned rotating heat pipe was studied by Schafer 
[Ref. 6]. The two-dimensional model was studied by Corley [Ref. 7], The 
same assumptions and boundary conditions, similar to those vised in the 
Nusslet analysis of film condensation on a flat wall, and based upon the 
analysis of Ballback [Ref. 2] were used for both. The more important of 
these assumptions are: 

1. steady state operation, 

2. film condensation, as opposed to drop wise condensation, 

3. laminar flow of the condensate film along both the fin and the 
trough, 

4. static balance of forces within the condensate, 

5. one-dimensional conduction heat transfer through the film thickness 
(no convective heat transfer in the condensate film), 

6. no liquid-vapor interfacial shear forces, 

7. no condensate subcooling, 

8. zero heat flux boundary conditions on both sides of the wall section 
(symmetry conditions), as shown in figure 5, 

9. saturation temperature at the fin apex, 

10. zero film thickness at the fin apex, and 

11. negligible curvature of the condenser wall. 
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Figure 3 shows the linear triangular finite element model developed 
by Purnomo [Ref. 1] for use in obtaining a two-dimensional Finite Element 
solution. 

The assumption that was used by Corley [Ref. 7] that the fin apex was 
at the saturation temperature of the working fluid was modified by 
Purnomo [Ref. 1], The value of the temperature at the apex of the fin was 
allowed to float and a parabolic temperature distribution was assumed along 
the fin surface. 

Purnomo’ s problem statement for the formulation of the Finite Element 
Method as shown in figure 5 is: 



jF + j£r_ 0 

dX 2 + dy 1 

with the following -boundary conditions: 



( 6 ) 



1 . along boundary 1 , -K dT/dn-h^T-Tsat ) 



2. along boundary 3, — K dTIdn-fyT-T*) 



dT 

3. along boundaries 2 and 4, -0 

dr i 



A detailed description of the numerical formulation is presented in his 
thesis. 

Davis used Constrained Function Minimization (CONMIN) as an 
optimization program. CONMIN is a FORTRAN program in subroutine form. 
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b) - !< ~ s (T - T a ) Along Boundary [3] 

c) !-i = 0 Along Boundaries C2] and C^] 



Fi9 «.'vf e * 5 ‘ 1 Di ^ er ® ntial Equation and Boundary Conditions Considered 
m the Analysis of Purnomo [Ref. 1]. 
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Vanderplaats [Ref. 11] developed the Control Program for Engineering 
Synthesis (COPES) as a main program to simplify the use of CONMIN. Davis’ 
computer program was written in subroutine form with SUBROUTINE ANALIZ 
(ICALC) as the main routine. The name ANALIZ is compatible with the COPES 
program and ICALC is a calculation control. Subroutine ANALIZ calls other 
subroutines as needed: 

1. the routine ’’COORD" used to define positions of system coordinate 
points 

2. the routine "FORMAF" used to formulate the Finite Element Method 
equations, 

3. the routine "BANDEC" as an equation solver for a symmetric matrix 
which has been transformed into banded form, and 

4. the routine "DPLORT" used to compute the roots of a real polynomial 
using a Newton -Rap hson derivative technique. 

Bo THE COMPUTER PROGRAM DEVELOPMENT 

The basis for the present analysis code is Davis’ [Ref. 9] two- 
dimensioned finite element program. Davis’ code was checked for validity as 
the first task undertaken in the development of this thesis. An error was 
discovered in calculating the fin condensate film thickness (AZS). In the 
initial cedculation of HDEN only, the cubed term was merely multiplied by 
three. The correct form of the equation is shown below: 

HDEN - -<J,z 3 /3 -6,z J /2 +z(T sat - 7\) 

The effect of this error was minimal since subsequent equations were 
correct, a 0.00016% difference in the condensate level was noted. 
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The next task undertaken was to adapt the analysis code to permit 
automated design analysis using ADS. Many modifications were made, some 
of which are mentioned here. The program was rewritten to include a main 
program from which ADS is called and subroutines to perform various 
mathematical functions. Subroutine ANALIZ was deleted since the double 
precision version of ADS (DADS) was used. Modifications were also made to 
generalize the code and minimize the changes needed when varying the 
number of finite elements used. 

A listing of the revised computer program is included as the 
Appendix. 

C. DESIGN OPTIMIZATION 

There are thirteen parameters that can be used as design variables. 
There are geometric or functional parameters of the rotating heat pipe or 
the properties of the working fluid or environment. The possible design 
variables, possible constraint functions, and the objective function are 
listed below in Fortran. This code can pursue a wide variety of design 
problems. 

The addition of fins by the designer increases the surface area which 
increases the heat transfer rate through the condenser wall. However, the 
addition of fins decreases the cross-sectional area through each fin for 
conduction and decreases the trough width which increases the condensate 
thickness in the trough. The increased condensate thickness decreases the 
heat transfer and if increased to the point of covering the fins it could 
dramatically reduce the heat transfer through the fin. 
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TABLE H. DESIGN VARIABLES 



DESIGN VARIABLES 



BFIN (fill height) 

CANGL (cone half angle) 

CLI (condenser length 
FANGL (fin half angle) 

HINF (external convective heat transfer coefficient) 
NFIN (number of fins) 

R2I (intermediate radius) 

RBASEI (inside radius of condenser base) 

RMP (rotational speed of the heat pipe) 

THICKI (condenser wall thickness) 

TINF (external temperature) 

TSS (saturation temperature) 

TZ (nodal point temperature) 



CONSTRAINT FUNCTIONS 



BOA (ratio of fin height over fin width) 
ZOA (ratio of trough width over fin width) 
DEL(NI) (condensate thickness) 



OBJECTIVE FUNCTION 



OBJ 



The purpose of the design study was to determine the fin height, 
number of fins, and fin half angle which would maximize the heat transfer 
rate. It was then decided that the design variables would be BFIN, FANGL, 
and NFIN. The number of fins was chosen vice the ratio of trough width 
to fin width (ZOA)' as was previously used. This decision was made since 
it is easier to' think in terms of the number of fins vice a ratio. Other 
potential design variables were held constant. The objective function to be 
maximized was OBJ=QT+QTF, the heat transfer through the fin plus the heat, 
transfer through the trough. 

The code was run with the three design variables using the internal 
scalar default parameters in ADS. The objective function was calculated 
using the input values of the design variables. The ADS program then 
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changed the first design variable keeping the other two design variables 
at the input values. The calculations were made aging yielding a different 
value for the objective function. The new value would then be compared 
to the previous value of the objective function. If the difference in the 
objective function was not greater then the internal scalar default value, 
the first design variable would be returned to its original value and the 
process repeated with the second design variable. If the difference in the 
objective function wag still not areater then the internal scalar default 
value, the second design variable would be returned to its initial value and 
the process repeated with the third design variable. 

Each time the program was run, the optimization code would choose 
BFIN as the design variable to change first as it .had the greater influence 
on the objective function. The remaining two variables would then either 
be kept constant or the number of fins would be maximized and the fin 
half angle minimized. Consistent results were not obtained with this method. 

To improve the results, the internal scalar parameters were modified. 
These modifications included the constraint tolerances, the absolute and 
relative convergence criteria, the absolute and relative change in the 
design variables, the absolute and relative change in the objective 
function, the minimum absolute value of the finite difference step when 
calculation gradients, and the initial relative move limit. Better results were 
obtained as seen in the higher value for the objective function. However, 
the results were still not consistent depending on the initial values used 
for the three design variables. At this point, it was decided to concentrate 
on one design variable and on the basis of the previous calculations the 
design variable chosen was BFIN. 
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The external surface temperature was set equal to the working fluid 
saturation temperature and the theoretical upper limit on the heat transfer 
was calculated for comparison. An assumption is made that there is no 

Q 

thermal resistance across the condensate, or the condenser wall. The upper 
limit of the heat transfer rate was predicted to be 69,492 BTTJ/HR using the 
following formula: 

- InrKT^TJ tt) 

where 

h = outside convective heat transfer coefficient (5000 BTU/HR'FT F) 
r = average outside radius of condenser wall (0.07373 FT) 

1 = condenser length (0.75 FT) 

T wa ji = temperature of the outside wall (100*F) 

T ffl = ambient temperature (60 *F) 

Based on engineering judgement certain constraints were placed on the 
design. These constraints were applied to the number of fins (not to 
exceed 400) and the minimum fin half angle (not to be less than 10 
degrees). The values were based on structural and manufacturing 
considerations. 
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IV. RESULTS 



A. INTRODUCTION 

The purpose of the design optimization was to- maximize the heat 
transfer rate. This was accomplished by using the computer code in 
conjunction with ADS. Numerical results are discussed below. 

B. CONSTRAINED OPTIMIZATION 

In the design problem undertaken to determine the optimum internal 
geometry for the maximum heat transfer, numerous runs were made for a 
condenser made of copper. This material has a thermal conductivity of 230 
BTU/HR’FT'F. The working fluid was water. 

Since the fin height (BFIN) was the design variable, the initial runs 
investigated whether there was an optimum fin height. Initially the fin half 
angle was held constant at 10 degrees and the number of fins was varied 
from 150 to 400 . In each case, for the optimum design, the fin height was 
maximized and the trough was eliminated, as seen in figures 6 and 7. 

The number of fins was then held constant and the fin half angle was 
varied from 10 to 25 degrees with the fin height remaining the design 
variable. Once again, the greatest heat transfer rate was achieved with the 
highest fin height for each number of fins (figure 8.) 

As seen in figure 9, the highest heat transfer rate achieved was for 
400 fins with a 10 degree fin half angle and a fin height of 0.0345 inches. 
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HEAT TRANSFER RATE ( Btu/»ir) 
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Figure 6. Heat Transfer Rate vs. Fin Height. 
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HEAT TRANSFER RATE ( btu/hr ) 




LEGEND 

□ = NUMBER OF RNS=500 
o = NUMBER OF FJNS=550 
a = NUMBER OF FiNS=4C0 



Figure 7. Heat Transfer Rate vs. Fin Height. 
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HEAT TRANSFER RATE (BTU/HR) 
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Figure 8. Heat Transfer Rate vs. Fin Half Angle (Optimum Fin Height). 



28 



63000 



62000 ” 
61000 * 
60000 - 
• 59000 - 
l 58000 '- 
57000 - 
56000 - 
55000 - 
54000 - 
53000 - 
52000 - 

siooo-L 




50 100 150 200 250 

NUMBER CF FINS 



300 



350 



400 



Figure 9. Heat Transfer Rate vs. Number of Fins (Optimum Design). 
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Figure 10 shows a plot of heat transfer rate versus fin half angle for 
a condenser with between 100 and 400 fins. The heat transfer rate, as a 
function of fin half angle for a constant fin height, increases as the fin 
half angle increases. Davis [Ref. 9] concluded that the heat transfer rate 
increased with an increase in fin half angle. This is correct if the fin 
height is kept constant. As the fin half angle increases, the surface area 
of the fin increases. The added surface area also has a thinner film of 
condensate on it which offers lower thermal resistance. The trough area 
decreases and the condensate film in the trough thickens, increasing the 
thermal resistance. This degradation does not offset the gain in the heat 
transfer rate caused by the fin. 

For external heat transfer coefficients varying from 1000-50,000 
BTU/HR»FT «F, the same optimum design geometry for a maximum heat 
transfer rate was obtained, which is stated in table III below. Figure 11 
shows the strong influence the external heat transfer coefficient has on 
the heat transfer rate. 

In figure 12, the effect of the rotating speed on the heat transfer 
rate is seen. As the rotational speed increases, the heat transfer rate 
increases, this is caused by an increase in the element heat transfer 
coefficient and a decrease in the condensate thickness on the fin which 
lowers the thermal resistance. 
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HEAT TRANSFER RATE (BTU/HR) 
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Figure 10. Heat Transfer Rate vs. Fin Half Angle (Constant Fin 
Height). 
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HEAT TRANSFER RATE (BTU/HR) 
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Figure 11. Heat Transfer Rate vs. Heat Transfer Coefficient. 
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HEAT TRANSFER RATE (btu/hr) 




LEGEND 

□ = RPM= 3600 
o = RPM= 3000 
a = RPM= 2500 



Figure 12. Heat Transfer Rate vs. Number of Fins (RPM Variation). 
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CONDENSATE THICKNESS (FEET) 




LEGEND 

□ = NUMBER OF FINS=100 
o = NUMBER OF FINS=200 
a = NUMBER OF nNS=3C0 
* = NUMBER OF F1NS=4-00 



Figure 13. Condensate Level vs. Position (100-400 Fins). 
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When the number of fins was increased from 100 to 400, maintaining 
the same fin half angle, the condensate level decreased with the increased 
heat transfer rate (figure 13). This occurred because the thinner film over 
the fins decreased the resistance across the film which raised the 
temperatures along the fin, which in combination with the lower height fins 
increased the temperature on the outside of the pipe. This increase in 
temperature brings the operation closer to the condensing limit. When the 
fin half angle was increased for a specified number and height of fins, the 
condensate level decreased due to the increased trough width (figure 14). 



TABLE m. OPTIMIZATION RESULTS FOR 
VARIOUS HEAT TRANSFER COEFFICIENTS 



OPTIMUM DESIGN 


Fin Height 


0.0345 inches 


Fin Half Angle 


10.0 degrees 


Number of Fins 


400 


^external 


HEAT TRANSFER RATE 


(BTU/HR*FT 2 *F) 


(BTU/HR) 


1000 


13765 j 


5000 


62252 


50000 


261003 



Figures 15 and 16 show the effect the increase of surface area has 
on the heat transfer rate. In figure 15, the fin height for 400 fins with a 
10 degree half angle is plotted against the heat transfer rate. As the fin 
height is increased up to a maximum value of 0.0345 inches the resulting 
design is a sawtooth. Figure 16 shows the effect adding more fins has on 
the heat transfer rate for a constant height fin. The greatest increase in 
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the heat transfer rate is seen from the addition of 100 fins from a smooth 



tube. 



In figure 17, the ratio of the actual surface area over the surface 
area for a smooth tube is plotted versus the number of fins. The fin 
height and the fin half angle cure both held constant. As expected, the 
surface area ratio increases in a linear manner as the number of fins 
increases. The ratio of the heat transfer rate over the heat transfer rate 
for a smooth tube is plotted for two different heat transfer coefficients, 
h=1000 BTU/HR»FT 2 *F and 5000 BTU/HR*FT 2 *F. An increase in the number 
of fins results in not only an increase in the area ratio but also an 
increase in the heat transfer ratio. The increase in the heat transfer ratio 
is greatest when going from a smooth tube to a tube with 100 fins. The 
heat transfer ratio increase is greater for the heat transfer coefficient 
equal to 5000 BTU/HR»FT *F. This is because the heat transfer coefficient 
has a direct effect on the heat transfer rate, that is. 



Both curves approach an asymptotic value. However, the curve with the 
lower heat transfer coefficient approaches this asymptotic value with a 
fewer number of fins. Additionally, in view of the relatively small increase 
in the heat transfer ratio by the addition of fins for the lower heat 
transfer coefficient, consideration should be given to the cost of 
manufacturing the fins versus the benefit derived by their addition. 
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CONDENSATE THICKNESS (FEET) 




LEGEND 

a = FIN ANGLE = 10 
o = FIN ANGLE = 15 
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+ = FIN ANGLE = 25 



Figure 14. Condensate Level vs. Position (400 Fins, 10-25 Degree Fin 
Half Angle). 
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HEAT TRANSFER RATE ( btu/hb ) 




Figure 15. Heat Transfer Rate vs. Fin Height (400 Fins). 
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HEAT TRANSFER RATE (BTU/HR) 
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Figure 16. Heat Transfer Rate vs. Number of Fins (Constant Fin 
Height). 
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RATIO 




LEGEND 
□ = H=5000 
o = H=10Q0 



Figure 17. Heat Transfer and Area Ratios vs. Number of Fins. 
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C. AUTOMATED DESIGN SYNTHESIS (ADS) 

This optimization project was done on the IBM mainframe using the 
ADS optimization program. As stated previously, ADS is a general purpose 
numerical optimization program with a variety of algorithms that can be 
used to tailor the solution. The solution is separated into three levels in 
ADS: strategy, optimizer, and one-dimensional search. For this problem the 
following combination of algorithms were used: 

1. Strategy: Sequential Linear Programming 

2. Optimizer: Modified Method of Feasible Directions for constrained 
minimization 

3. One-Dimensional Search: Golden Section Method followed by 

polynomial interpolation 

The strategy used linearizes a nonlinear problem by a first order 
Taylor series expansion of the objective and constraint functions. The 
solution to this linear approximation is obtained. The problem is linearized 
again about this point and the new problem is solved with the process 
being repeated until a precise solution is achieved. 

The optimizer chosen is used to find a search direction which will 
minimize the objective function while maintaining feasibility. 

The combination of strategy, optimizer, and one-dimensional search 
chosen is not the only one available, nor is it necessarily the most 
efficient. It did yield results that were maximized and were within the 
constraint tolerances. 

The ADS optimization program is complicated by the numerous internal 
parameters which must be changed to obtain an optimal design. 



41 



Complications arise when there is a vast difference in the scales of the 
design variables. The design variables themselves must be scaled which is 
further complicated when the variable itself covers a wide range. ADS also, 
in certain cases, allows for constraints to be violated. In some instances 
this might be acceptable but not in this case. 

ADS does not have a scoping mechanism, that is the ability to 
decrease the rate of change of the design variable, and therefore to obtain 
a precise answer the internal parameters must be changed repeatedly. ADS 
also does not recognize integers, all numbers are real therefore, depending 
on the answer given, it may be necessary to round up or down. 
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V. CONCLUSIONS 



1. For an independent increase in fin half angle, rotational speed or 
number of fins an increase is seen in the heat transfer rate. As the 
parameters increase, the heat transfer rate levels off at the theoretical 
maximum heat transfer rate for the heat pipe. A decrease in the fin half 
angle with a corresponding increase in the fin height increases the heat 
transfer rate. If the fin half angle is decreased while the fin height is 
kept constant, then the heat transfer rate decreases. 

2. Maximum heat transfer occurs for the same fin geometry 
regardless of the external heat transfer coefficient. For a specific 
condenser radius, as many fins as possible should be machined with a 
minimum fin half angle at the maximum fin height. 

3. The computer code can be used for single analysis or the 
automated design of an internally finned rotating heat pipe. 

4. The benefit of adding fins is dependent on the external heat 
transfer coefficient. Consideration of the cost of manufacturing the fins 
versus the increase in the heat transfer rate should be made. 
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VI. RECOMMENDATIONS 



1. Analyze different shaped fins including rectangular and curved, 

2. Modify the code to allow for silmutaneous variations of more than 
one variable. 

3. Use different working fluids and heat pipe materials to see if a 
different internal geometry occurs for the maximum heat transfer rate. 

4. Modify the code to use the DOT optimization program vice ADS. 
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APPENDIX: PROGRAM LISTING 



PROGRAM OLSON 

******************************************************** 



* * 

* ANALYSIS OF ROTATING HEAT PIPE , USING TRIANGULAR * 

* ELEMENT MODEL * 

* COMPILED BY MAJOR IGNATIUS . S . PURNOMO IN JUNE 1978 * 

* * 

* MODIFIED TO PERMIT NUMERICAL OPTIMIZATION * 

* USING COPES/CONMIN * 

* BY LCDR WILLIAM A. DAVIS, JR. IN SEPTEMBER 1980 * 

* MODIFIED TO PERMIT USE OF THE ADS OPTIMIZATION * 

* ROUTINE BY LT. G.L. OLSON IN JUNE 1992 * 

* 

* * 



******************************************************** 



CHARACTER* 20 NAME 

. COMMON/AD S/D F( 21) ,G(10) ,IDG(100) , IGRAD, INFO, IOPT, IONED, IPRINT, 
rISTRAT, IWK( 2000) , IZ( 30) ,OBJ,S( 2 ) , VLB ( 2) ,VUB( 2) ,W( 21 , 30) ,WK( 5000 ) , 
: NCOLA , NCON , NDV , NGT , NRA , NRIWK , NRWK 



COMMON/OLL I E/A ( 200,50) , AMTOT( 200),APS,B(3) , BFIN , BOA , BVIN , C ( 3 ) , 

: CANGL , CF ( 200 ) , CK , CLI , COF ( 5 ) ,CW( 200) , DEL ( 200) , DMDOT ( 200),EA(3,3), 
: EPS ( 200) , EZERO , F ( 200,1) , FANGL, FNOBJ( 100 ) ,H( 200 ) , HINF , HZ ( 200 ) , 

: QB ( 200 ) , QINC ( 200 ) ,QTINC(200) , QTOT , QTOTAL (100) ,R(200),RB(200), 
iRBASEI , R2l , RHOF( 200 ) ,ROOTI(4) , RPM , ROOTR ( 4 ) ,T(200) , TALFA , TB ( 20 0 ) , 
:TCC( 200 ) , TE ( 200 ) , THICK , THICKI , TIB ( 200 ) , TINF , TS ( 200 ) ,TSAT,TSS, 
:TT( 200 ) ,TZ,UF( 200 ) ,X( 200) ,XCOF(5) ,XPLOT(200) ,Y(200),Z(200) ,ZOA, 

: DOBF , DOTH , I COR ( 200,3) , I FF , JINT , JLC , JTC , KFF ( 50) , KFIN ( 50 ) , KT , NBAN , 
: NEL , NFIN , NSNP 



GUIDE TO FORTRAN VARIABLE NAMES 



AFOVAS 


FIN AREA/SMOOTH AREA 




ALFA 


FIN HALF ANGLE (RADIANS) 




BFIN 


HEIGHT OF FIN (INCHES) 




BOA 


TANGENT OF THE FIN HALF ANGLE 




BVIN 


HEIGHT OF FIN (FEET) 




CAL FA 


COSINE OF ALFA 




CANGL 


CONE HALF ANGLE (DEGREES) 




CBASE 


INSIDE CIRCUMFERENCE OF CONDENSER 


( FEET) 


CEXIT 


INSIDE CIRCUMFERENCE AT CONDENSER 


EXIT (FEET) 


CF 


THERMAL CONDUCTIVITY OF CONDENSATE 


FILM ( BTU/HR FT F 


CL 


CONDENSER LENGTH (FEET) 




CLI 


CONDENSER LENGTH (INCHES) 




CPHI 


COSINE OF PHI 




CRIT 


CONVERGENCE CRITERION 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



DEL FILM THICKNESS 

DF GRADIENT OF OBJECTIVE 

DIV FLOATING POINT VALUE OF NDIV 

DMTOT CONDENSATE MASS FLOW RATE 

EPS TROUGH WIDTH INCLUDING INCREMENTAL CHANGE 

EPSEX TROUGH WIDTH AT CONDENSER EXIT 

EPSO TROUGH WIDTH AT START OF CONDENSER 

EZERO FIN BASE WIDTH 

F FORCE VECTOR OF SYSTEM 

FANGL 'FIN HALF ANGLE (DEGREES) 

G CONSTRAINT VALUES ASSOCIATED WITH CURRENT DESIGN 

H CONVECTIVE HEAT TRANSFER COEFFICIENT ( BTU/HR FT2 F) 

HFG LATENT HEAT OF VAPORIZATION ( BTU/LBM ) 

IDG CONSTRAINT TYPE IDENTIFIER 

IEL THE ELEMENT NUMBER 

IFF NO. OF ROWS MINUS ONE OF THE UPPER TRIANGULAR FIN 

IFIN EQUALS 0 FOR COPPER, AND 1 FOR STAINLESS STEEL 

IFLUID EQUALS 0 FOR WATER, AND 1 FOR FREON 

I GRAD GRADIENT CALCULATION IDENTIFIER 

INFO CONTROL PARAMETER 

IONED ONE DIMENSIONAL SEARCH IDENTIFIER 

IOPT OPTIMIZER IDENTIFIER 

IPRINT A FOUR DIGIT PRINT CONTROL 

ISTRAT OPTIMIZATION STRATEGY IDENTIFIER 

JINT NO. OF COLUMNS PLUS ONE BELOW TRIANGULAR FIN 
JLC NUMBER OF SYSTEM NODAL POINT LOCATED AT 

THE CENTER OF SYSTEM COORDINATE 
JTC NUMBER OF SYSTEM NODAL POINT LOCATED AT 

THE JUNCTION OF THE SYMMETRY BOUNDARY AND 
THE LINE OF INTERSECTION BETWEEN THE FIN 
AND THE CONDENSER WALL 

KFF NUMBER OF SYSTEM NODAL POINTS LOCATED ALONG 

THE FIN CONVECTIVE BOUNDARY 
KFIN NUMBER OF SYSTEM NODAL POINTS LOCATED ON THE 
SYSMMETRIC BOUNDARY OF TRIANGULAR FIN SECTION 
NOT COUNTING POINTS AT BASE AND APEX 
KT NUMBER OF COLUMNS WITHIN THE TROUGH WALL SECTION 

NBAN SYSTEM BAND WIDTH 

NBOTF LAST ELEMENT AT BOTTOM SIDE 

NBOTI FIRST ELEMENT AT BOTTOM SIDE 

NCON NUMBER OF CONSTRAINTS 

NDOBF NUMBER OF ROWS WITHIN THE FIN 

NDOTH NUMBER OF ROWS WITHIN THE TROUGH 
NDIV NUMBER OF INCREMENT 

NDV NUMBER OF DESIGN VARIABLES 

NEFB ELEMENT NUMBER AT BASE OF FIN 

NEL NUMBER OF ELEMENTS 

NEST ELEMENT NUMBER AT END OF TROUGH 
NRA NUMBER OF ROWS IN ARRAY A 

NRWK DIMENSIONAL SIZE OF WK 

NSNP NUMBER OF SYSTEM NODAL POINTS 

OBJ VALUE OF THE OBJECTIVE FUNCTION ASSOCIATED WITH X 

OMEGA ROTATIONAL SPEED OF HEAT PIPE (RAD/HR) 

PHI CONE HALF ANGLE (RADIANS) 

PI PI 

R2I DISTANCE FROM CENTERLINE OF THE HEAT PIPE TO HALF THE F 

HEIGHT 

RBASE INSIDE RADIUS OF CONDENSER BASE (FEET) 

RBASEI INSIDE RADIUS OF CONDENSER BASE (INCHES) 

REXIT INSIDE RADIUS OF CONDENSER EXIT (FEET) 
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RPM REVOLUTIONS PER MINUTE 

S VECTOR OF DESIGN VARIABLES 

SALFA SINE OF ALFA 

SPHI SINE OF PHI 

SURFAR SURFACE AREA 

THICK CONDENSER WALL THICKNESS (FEET) 

THICKI CONDENSER WALL THICKNESS (INCHES) 

TPHI TANGENT OF PHI 

TZ AMBIENT TEMPERATURE 

UF VISCOSITY 

VLB LOWER BOUNDS ON THE DESIGN VARIABLES ’ 

VUB UPPER BOUNDS ON THE DESIGN VARIABLES 

WK REAL WORK ARRAY 

ZFIN NUMBER OF FINS 

ZOA RATIO OF TROUGH WIDTH TO FIN BASE WIDTH 

ZSTAR SURFACE LENGTH OF THE FIN MINUS THE SURFACE LENGTH 
COVERED BY THE CONDENSATE IN THE TROUGH 
ZZERO SURFACE LENGTH OF FIN 



PRINT*, 'INPUT FILE NAME' 
READ*, NAME 
OPEN ( 10 , FILE-NAME ) 
OPEN(15,FILE-'/HTPIPE OUTPUT') 
OPEN( 14, FILE- '/DUMP OUTPUT') 
OPEN( 13, FILE-'/GRAPH OUTPUT') 



THE FOLLOWING READS INPUT DATA, PERFORMS HEAT 
TRANSFER ANALYSIS, AND PRINTS RESULTS. 



***** INPUT MODE ***** 



ELEMENT CONNECTIVITIES 

READ (10,420) NEL , NSNP , NBAN , I FLUID , I FIN 
WRITE (15,430) NEL, NSNP, NBAN 
WRITE ( 15,435) I FLUID , I FIN 
WRITE (15,436) 

WRITE (15,437) 

READ (10,440) ( I CL , ( ICOR( IEL,I) ,1-1,3) ,IEL=1, NEL ) 
WRITE (15,450) 

THE CONDENSER GEOMETRY 

READ (10,460) CLI , CANGL , RBASEI , R2I , THICKI ,BFIN,TZ 
WRITE (15,470) CLI , CANGL , RBASEI , R2l , THICKI , BFIN , TZ 
READ (10,480) NDIV , NEST , NEFB , NBOTI , NBOTF 
WRITE (15,490) NDIV, NEST, NEFB, NBOTI , NBOTF 

DATA FOR RUNNING 

READ (10,500) RPM, TSS , TINF , HINF 
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no * * * * * * * *o**o 



WRITE (15,510) RPM, TSS , TINF , HINF 

THE CONVERGENCE CRITERIAN 

READ (10,520) CRIT 
WRITE (15,530) CRIT 

INTERNAL FIN GEOMETRY 

READ (10,540) FANGL , NFIN 
WRITE (15,550) FANGL, NFIN 
WRITE ( * , * ) FANGL, NFIN 
READ (10,560) IFF 
WRITE (15,570) IFF 

READ (10,580) (KFIN(I) ,KFF(I) ,1*1, IFF) 

READ (10,590) NDOBF , NDOTH , JTC , JLC , JINT , KT 
NHB-NEFB/2 
NBF-NBOTF+1 
DOBF - FLOAT (NDOBF) 

DOTH * FLOAT (NDOTH) 

WRITE (15,600) I COR ( NBOTI , 2 ) , I COR ( NEFB , 1 ) , I COR ( NEST , 1 ) , 

1 I COR ( NBOTF , 1 ) 

SET CONSTRAINTS 
NRA-21 
NCOLA-30 
NRWK- 5000 
NRIWK-2000 
NDV- 1 
NCON- 2 
I GRAD- 0 

INITIAL DESIGN 
S( 1 )- BFIN 

BOUNDS 

VLB ( 1 )- 0.0000001 
VUB ( 1 ) - 0 . 7 5 

IDENTIFY CONSTRAINTS 
NONLINEAR CONSTRAINT 
IDG ( 1 )-l 
C LINEAR CONSTRAINT 

IDG ( 2 ) -2 

* 

PRINT*, 'INPUT THE VALUES FOR ISTRAT, IOPT, IONED AND IPRINT' 
READ*, ISTRAT, IOPT, IONED, IPRINT 



C INITIALIZE COUNTER 

NO-O.O 

C CHANGE THE INTERNAL PARAMETERS 

INFO— 2 

CALL DADS( INFO, ISTRAT, IOPT, IONED, IPRINT, IGRAD, NDV, NCON, S, VLB, 
: VUB , OBJ , G , IDG , NGT , IZ , DF , W , NRA, NCOLA, WK , NRWK , IWK , NRIWK ) 

IWK ( 2 ) =0 
IWK ( 3 ) — 2 0 0 
IWK ( 5 ) =4 
IWK ( 7 ) =500 
WK(3)— 0.5 
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WK(6)— 0.01 

WK(8)-0.05 

WK( 9 )-0 . 50 

WK(10)-0.05 

WK(11)=0.005 

WK( 13 )=0 . 001 

WK (14J-0.0001 

WK ( 21 ) =0 . 004 

WK( 22 )-0 . 002 

WK( 26 )=0 .004 

WK(37)=0. 0000000001 

10 CALL DADS ( INFO , I STRAT , I OPT , IONED , I PRINT , I GRAD , NDV , NCON , S , VLB , 
: VUB , OBJ , G , IDG , NGT, IZ , DF , W, NRA, NCOLA, WK , NRWK , IWK , NRIWK ) 

IF (INFO.EQ.0) GO TO 360 

***** EXECUTION MODE ***** 



NO-NO+1 

CONVERT UNITS OF ALL DIMENSIONAL PARAMETERS 
TO FEET. CONVERT UNITS OF ANGLES TO RADIANS. 
R2l-RBASEI-( S( 1 )/2 ) 

CL-CLI/12.0 
R2-R2I/12.0 
RBASE-RBASEI/12 . 0 
BVIN-S ( 1 )/12 . 0 
DIV-FLOAT(NDIV) 

PI-3. 14159265358979 
PHI-2. 0*CANGL*PI/360.0 
SPHI-SIN ( PHI ) 

CPHI-COS( PHI ) 

TPHI-TAN( PHI ) 

DELX-CL/DIV 
CBASE-2 . 0*PI*RBASE 
REXIT=RBASE+CL*TPHI 
CEXIT-2 . 0*PI *REXIT 
THICK-THICKI/12 . 0 
ALFA-FANGL*2 . 0*PI/360 . 0 
SAL FA-SI N ( AL FA ) 

CALFA-COS ( ALFA ) 

TAL FA- TAN (ALFA) 

EZERO-2 . 0* ( S ( 1 )/12 . 0 ) *TALFA 

ZOA— ( ( CBASE- ( EZERO*NFIN ) ) /NFIN ) /EZERO 

BOUNDARY CONDITIONS AND TEMPERATURE ESTIMATES 
ALONG THE FIN BOUNDARY 

DO 20 NTINF-NBOTI ,NBOTF 
20 TS ( NTINF ) -TINF 
DO 30 NNT-NBF , NEL 
TS ( NNT ) =0 . 0 
30 H ( NNT ) -0 . 0 

DO 40 IGT-1 , NEST 
IE-ICOR( IGT , 2 ) 

40 T ( I E ) =TZ 

IG— ICOR( NEST, 1 ) 

T( IG ) — TZ 

OMEGA IS IN RADIANS/HOUR 
OMEGA=RPM*2 . 0*PI*60 . 0 
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DO 50 KL-NBOTI , NBOTF 
50 H{KL)-HINF 
HIFN-HINF 
TSAT-TSS 
EPSO=ZOA*EZERO 
BOA=* TAL FA 

SURFAR-NFIN* ( 2 . 0* ( S ( 1 )/( 12*CALFA) )+EPSO) 

EPS EX— ( CEXIT— ( NFIN*EZERO ) )/NFIN 

BETA- ( EPSEX-EPSO )/DIV 

ZZERO— ( S ( 1 )/12 )/CALFA 

AFOVAS- ( ZOA+ ( 1 ./SALFA) )/( 1 .+ZOA) 

ZA-0 c 0 

DO 60 NSAT-1 , NEST 
60 TS ( NSAT ) -TSAT 

TSOLID- ( TSAT+TINF )/2 . 0 
C TEMPORARY CHANGE - TFILM 

ASMOOTH-O . 0 
ACASE-0. 0 
QT-0.0 
OBJ-O.O 
QTl-0.0 
QTF-0.0 
QTRF-0.0 
QTOT-O . 0 
DMTOT-O . 0 
NK-NDIV+1 
DO 350 NI-1 , NK 

C R IS THE INCREMENTAL CHANGE IN THE RADIUS OF THE CONDENSER 

R( NI ) -R2+NI *DELX*SPHI 
RB ( NI ) -RBASE+NI *DELX*SPHI 

C EPS IS THE INCREMENTAL CHANGE IN THE TROUGH WIDTH 

EPS ( NI ) -EPSO+NI *BETA 
APS-EPS ( NI ) 

NODAL POINT COORDINATES 

CALL COORD 
65 Z { 1 ) -ZA 

DO 70 I ZEL-1 , NEFB 
NA-ICOR( IZEL , 1 ) 

NB- 1 COR ( I Z EL , 2 ) 

XE— X ( NA ) -X ( NB ) 

YE= Y ( NA ) -Y ( NB ) 

ELZ-SQRT ( XE* *2+ YE* * 2 ) 

70 Z( IZEL+1 )=Z( IZEL)+ELZ 

XZB-X( I COR ( NHB , 1 ) ) -X( ICOR( 1 , 2 ) ) 

YZB-Y ( ICOR( NHB , 1 ) )-Y( ICOR( 1,2 ) ) 

ZB-SQRT( XZB**2+YZB**2 ) 

ZC-ZZERO 
IM-1 

PARABOLIC TEMPERATURE DISTRIBUTION ALONG THE FIN 
BOUNDARY, USING LAGRANGE INTERPOLATION 

80 TPl-T ( I COR (1,2) ) 

TP2-T ( ICOR ( NHB , 1 ) ) 

TP3-T ( I COR ( NEFB , 1 ) ) 

APl=TPl/( ZB*ZC ) 

AP2-TP2/( ZB*( ZB-ZC) ) 

AP3=TP3/( ZC*( ZC-ZB) ) 
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BPl — (ZB+ZC)*APl 
BP2— ZC*AP2 
BP3— ZB*AP3 
A1-AP1+AP2+AP3 
B1-BP1+BP2+BP3 
TC-0 . 0 

DO 90 NY=1 , NEST 

90 TC-TC+T( ICOR(NY, 2) ) 

AY-FLOAT (NY+1 ) 

TF- ( TC+T ( I COR ( NY , 1 ) )+AY*TS(NY) )/(2.0*AY) 

SOLID-FLUID PROPERTIES 
WATER PROPERTIES 
I F( I FLUID . EQ . 1 ) GO TO 91 

HFG-1093.88-0.5703*TS(1)+0.00012819*(TS(1)**2) 

1-0.0000008824*(TS(1)**3) 

RHOF ( NI ) -6 2 .774-0 . 00255698*TF-0 . 000053572*TF* *2 
CF(NI )-0. 3034+0. 000738927*TF-0. 00000147 321 *TF* *2 
UF(NI )-0 . 001397-0 . 000014669*TF+0 . 0000 0006 312 53 *TF** 2-0 .00000 
100000976569* TF**3 

FREON PROPERTIES 

91 IF( I FLUID . EQ. 0 ) GO TO 92 

HFG-69. 5459-0. 0156011 *TS ( 1 ) -0 . 000455294* ( TS ( 1 )* *2 ) +0 . 00000104144* ( 
1TS( 1 ) **3 ) 

RHOF ( NI ) -102. 059-0. 025364 *TF-0. 000502649* (TF**2) +0. 0000 01 3 5407* (TF 
1**3) 

CF(NI ) -0.05948 58-0. 000429765*TF+0. 00000 348218*TF* *2-0. 00000001 04 16 
18*TF**3 

UF (NI )-0. 00078-0. 00000 525 *TF+0 ,0000000125*TF**2 

92 UF(NI )-3600*UF(NI ) 

IF( IFIN.EQ. 1 ) GO TO 93 

CW(NI ) -231 .7772-0 .022 22 *TSOLID 

93 IF( IFIN.EQ. 0 ) GO TO 94 

CW ( NI ) -8 . 776+0 . 00 265 *TS OLID 

94 CK-CW(NI) 

CONST-RHOF ( NI ) ** 2*OMEGA* * 2*HFG*CPHI *CALFA*R( NI ) 

INITIAL FILM THICKNESS 
IF (NI.GT.l) GO TO 100 

DEL(1)-1.107*( ( ( TSAT-TINF ) *CF(NI )/(UF(NI ) *HFG ) ) ** . 25) * ( (UF(NI )/( 
lRHOF ( NI ) * OMEGA ) ) * * 0 . 5 ) 

100 CONTINUE 

AVERAGE ELEMENT CONVECTIVE COEFFICIENT ALONG 
THE FIN BOUNDARY 

ZSTAR-ZZERO-DEL ( NI )/CALFA 
AZZ-DEL ( NI ) /SALFA 
ZZ-ZSTAR 

HDEN- ( ( -Al *ZZ * *3)/3.0-(Bl*ZZ**2)/2.0)+ZZ*( TSAT-T ( 1 ) ) 

AZS—ABS ( 4 *CF ( NI ) *UF(NI ) *HDEN/CONST ) **0 .25 
HAC-0.0 

DO 190 IEL-1 , NEFB 
AZ-Z ( IEL ) 

BZ—Z ( IEL+1 ) 

IF ( ZSTAR.LE.BZ ) GO TO 110 
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GO TO 120 

110 IF (HAC.NE.0.0) GO TO 180 
BZ-ZSTAR 

120 IF (IEL.NE.l) GO TO 130 
AK-( BZ-AZ )/5.0 
ZZ-AK 
GO TO 140 

130 AK-( BZ-AZ )/4 . 0 
ZZ-AZ 

140 ZEL»4*AK 

DO 150 NH-1,5 

HDEN= ( -1 . 0 * ( Al *ZZ * * 3/3 . 0+Bl*ZZ**2/2 . 0 ) )+ZZ* ( TSAT-T ( 1 ) ) 
HZ(NH)-ABS(CF(NI )**3*CONST/( 4 *UF ( NI ) *HDEN ) )**0.25 
150 ZZ-ZZ+AK 

C AVERAGE H USING SIMPSONS RULE 

CONH=-AK*(HZ( 1 ) + 4*HZ( 2 )+2*HZ( 3 ) + 4*HZ( 4 )+HZ( 5 ) )/( 3*ZEL) 
IF (ZSTAR.EQ.BZ) GO TO 160 
H( IEL)-CONH 
GO TO 190 
160 AZ-ZSTAR 

HAZ-CONH*(AZ-Z( IEL) ) 

DELA-AZS 
170 BZ-Z(IEL+1) 

DELB™ ( BZ-ZSTAR ) *AZZ/( ZZERO-ZSTAR) 

DELZ- ( DELA+DELB ) /2 . 0 
HAC*( BZ-AZ ) *CF( NI )/DELZ 
H( IEL ) ™ ( HAZ+HAC )/( BZ-Z ( IEL ) ) 

GO TO 190 
180 AZ-Z(IEL) 

DELA-DELB 
HAZ-0 .0 
GO TO 170 
190 CONTINUE 

NETI-NEFB+1 
DO 200 IEL-NETI ,NEST 
200 H( IEL)»CF(NI ) /DEL ( NI ) 

ENTRY INTO THE FINITE ELEMENT SOLUTION 

CALL FORMAF 

CALL BANDEC ( NSNP , NBAN , 1 ) 

THE TEMPERATURE DISTRIBUTION 

DO 210 NT=1 , NSNP 

210 T ( NT ) =F ( NT , 1 ) 

TIB ( NI ) *T ( ICOR( NBOTI , 2 ) ) 

TT ( NI ) =T( ICOR( NEFB , 1 ) ) 

TE ( NI ) =*T ( ICOR( NEST, 1 ) ) 

TB(NI )»T( I COR ( NBOTF , 1 ) ) 

TTS-0 .0 

DO 220 NS=1,NSNP 
220 TTS-TTS+T(NS) 

PN-FLOAT(NS ) 

TSOLID=TTS/PN 

Q AT THE BOTTOM SIDE 



QBI=0 . 0 

DO 230 IBEL=NBOTI , NBOTF 
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NKA-ICOR( IBEL,1) 

NKB-ICOR( IBEL, 2) 

XB-X(NRA)-X(NKB) 

YB-Y(NKA)-Y(NKB) 

ELB-SQRT ( XB* * 2+YB* * 2 ) 

QBI*QBI+(T(NKA)+T(NKB)-2*TS( IBEL) ) *ELB*H( IBEL)/2 
QB(NI )=QBI*DELX 



ITERATION UNTIL CONVERGENCE CRITERIA IS MET 



IF ( IM.EQ. 1 ) GO TO 240 
QJ-QBI 
GO TO 250 
240 QI-QBI 
IM-2 

GO TO 80 

250 AQ-ABS(QJ-QI )/QJ 

IF (AQ.LE.CRIT) GO TO 260 

QI-QJ 

GO TO 80 

260 DMDOT(NI )-2 . *QBI*DELX/HFG 
DMTOT=»DMTOT+DMDOT( NI ) 

Cl*RHOF ( NI ) * * 2 * OMEGA* * 2 *R ( NI ) *SPHI/( 3*UF(NI) ) 
XCOF(l)— DMTOT 



■ 0.0 
■ 0.0 

■C1*EPS ( NI ) 
■Cl*TALFA 



XCOF ( 2 ) ' 

XCOF ( 3 ) ' 

XCOF ( 4 ) ’ 

XCOF ( 5 ) 1 
M-4 

CALL DPOLRT ( M , IER ) 
IF (ROOTR(l) .GT.0.0) 
IF ( ROOTR ( 2 ) . GT .0.0) 
IF ( ROOTR ( 3 ) .GT. 0 . 0 ) 
IF ( ROOTR ( 4 ) .GT.0.0) 



GO 

GO 

GO 

GO 



TO 

TO 

TO 

TO 



270 

280 

290 

300 



THE CONDENSATE THICKNESS 



270 DEL ( NI+1 ) "ROOTR ( 1 ) 

GO TO 310 

280 DEL ( NI+1 )-ROOTR( 2) 

GO TO 310 

290 DEL ( NI+1 ) »ROOTR ( 3 ) 

GO TO 310 

•300 DEL ( NI + 1 ) =ROOTR ( 4 ) 

310 QEL-0.0 

IF (NI.NE.l) GO TO 320 

Q FROM THE TOP SIDE 



Q THROUGH FIN 



QEL-0.0 

320 DO 330 IQEL-1 , NEFB 
KA-ICOR( IQEL , 1 ) 

KB-ICOR( IQEL, 2 ) 

XQEL-X( KA)-X( KB) 

YQEL=Y ( KB ) -Y ( KA ) 

ELM-SQRT( XQEL**2+YQEL**2 ) 

QEL=»QEL+( 2 *TS ( IQEL ) -T ( KA) -T( KB ) ) * ELM*H ( IQEL )/2 
330 CONTINUE 
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QINC(NI )*QEL*DELX 
AMTOT(NI )-DMTOT 
QET«QEL*DELX*NFIN*2 
QT-QT+QET 

QA-QBI *DELX*NFIN*2 
Q THROUGH TROUGH 
QTRF-0 . 0 

DO 340 IQEL-NEFB+l ,NEST 
KA* I COR ( I QEL , 1 ) 

KB* I COR ( IQEL , 2 ) 

XQEL-X( KA ) -X ( KB ) 

YQEL-Y ( KB)-Y( KA) 

ELM-SQRT(XQEL**2+YQEL**2 ) 

QTRF-QTRF+( 2*TS ( IQEL ) -T( KA)-T( KB ) ) *ELM*H( IQEL )/2 . 0 
340 CONTINUE 

QTINC ( NI ) *QTRF*DELX 
QTOTAL ( NI ) -QINC ( NI ) +QTINC ( NI ) 

QTRFT=QTRF*DELX*NFIN*2 . 

QTF-QTF+QTRFT 
QTOT-QTOT+QA 

ASMOOTH-2*PI*RB(NI ) *DELX+ASMOOTH 

ACASE«ACASE+NFIN*( ( ( 2*S( 1 ) )/( 12*CALFA) )+EPSO) *DELX 
350 CONTINUE 

EVALUATE OBJECTIVE FUNCTIONS AND CONSTRAINTS 
OBJ— (QTOT) 

WRITE ( 15 , * ) 'OBJECTIVE FUNCTION * f , OBJ , ' BFIN * ' , S(l) 

FIRST CONSTRAINT IS TO ENSURE THE RATIO ZOA IS NOT NEGATIVE 
G( 1 )— ( ( ( CBASE- ( 2*NFIN*S ( 1 ) *TALFA/12 ) )/NFIN )/( S ( 1 ) *TALFA/6 ) ) 
G(1)-1000.0*G(1) 

THE SECOND CONSTRAINT IS TO ENSURE THE CONDENSATE LEVEL IS NO 
GREATER THAN THE FIN HEIGHT 
G( 2 ) — ( ( S ( 1 )/12 . 0 ) -DEL ( NI ) ) 

XPLOT ( NO ) *S ( 1 ) 

FNOBJ(NO)— OBJ 
ARATIO-ACASE/ASMOOTH 
WRITE ( 13 , * ) XPLOT(NO) ,FNOBJ(NO) 

WRITE (14,*) ARATIO, FNOBJ(NO) 

GO TO 10 
360 CONTINUE 
C 

BFIN=S(1) 

C ***** OUTPUT MODE ***** 

WRITE (15,630) 

DO 370 NR=1,NK 

370 WRITE (15,640) NR , QINC ( NR ), QTINC ( NR ), QTOTAL ( NR ) 

WRITE (15,650) QT , QTF 

WRITE (15,660) CL I , CANGL , RBASEI , R2 I , THICK I , BFIN , RPM , TSS , TINF , 
1RIT, FANGL , ZOA, IFF 
WRITE (15,661) AFOVAS 

WRITE (15,670) BOA, ZOA,NFIN, BVIN, SURFAR 
WRITE (15,680) 

DO 380 NP=1 , NSNP 

TCC( NP) = . 5555555* ( T( NP) -32 ) 

380 WRITE (15,690) NP , X ( NP ) , Y( NP ) , T( NP ) , TCC( NP ) 

WRITE (15,700) 

DO 390 KKL=*1 , NBOTF 
NKX-ICOR( KKL, 1 ) 

NKY=ICOR( KKL, 2 ) 
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XP-X(NKX)-X(NKY) 

YP-Y(NKX)-Y(NKY) 

EXY-SQRT(XP**2+YP**2) 

QEP-ABS ( ( T( NKX ) +T( NKY ) -2 *TS ( KKL ) ) *EXY*H ( KKL ) /2 . 0 ) 

QEP-QEP*DELX 

390 WRITE (15,710) KKL , H ( KKL ) , EXY , QEP 
WRITE (15,720) CRIT 

WRITE (15,730) HFG , NFIN , H ( NBOTF ) , TSAT , RPM , QTOT , QT , FANGL 
WRITE (15,734) 

WRITE (15,735) ROOTR(l) , ROOTI ( 1 ) 

WRITE( 15,735) ROOTR ( 2 ) , ROOTI ( 2 ) 

WRITE( 15,735) ROOTR ( 3 ) , ROOTI ( 3 ) 

WRITE( 15,735) ROOTR ( 4 ) , ROOTI ( 4 ) 

WRITE (15,740) 

DO 400 NR-1,NK 

400 WRITE (15,750) NR , DEL ( NR ) , QB ( NR ) , AMTOT ( NR ) , TIB ( NR ) , TT ( NR ) , TE ( NR ) , 
1TB (NR) 

WRITE (15,760) 

DO 410 NG-1 , NDIV , 2 

410 WRITE (15,770) NG, CW( NG ) , CF ( NG ) , RHOF ( NG ) , UF ( NG ) , EPS ( NG ) , R( NG ) , 
lQINC(NG) 

RETURN 



412 FORMAT ( 8X, E12 . 5 , 8X, E12 . 5 ) 

420 FORMAT (515) 

430 FORMAT ( /2X , 1 5HNO . OF . ELEMENTS- , 1 5 , 10X , 18HNO . OF . SYSTEM N.P.-, 

115 ,/, 2X, 13HNO .OF BANDED-, 15) 

435 FORMAT ( /2X, ' I FLUID-' , 15 , 10X, ' I FIN-' , 15 ) 

436 FORMAT (2X,'IFLUID - 0 FOR WATER, AND 1 FOR FREON') 

437 FORMAT (2X,'IFIN - 0 FOR COPPER, AND 1 FOR STAINLESS STEEL') 

440 FORMAT ( 4 1 5 ) 

450 FORMAT (/2X, 7HELEMENT, 10X, 3HNP1 , 14X, 3HNP2 , 15X, 3HNP3 ) 

460 FORMAT (7G10.5) 

470 FORMAT (4X,5HCLI- , E12 . 5 ,/, 4X, 7HCANGL- , E12 . 5 , / , 4X , 8HRBASEI . = , El 2 . 
15 ,/, 4X, 5HR2I- ,E12.5,/,4X, 8HTHICKI - , E12 . 5 , / , 4X , 6HBFIN= ,El2.5,/,4 
2X , 4HTZ- , El 2 . 5 ) 

480 FORMAT (5l5) 

490 FORMAT (4X,6HNDIV- , 110 ,/, 4X, 6HNEST- , I 1 0 , / , 4X , 6HNEFB- ,I10,/,4X,7 
1HNBOTI- , 110 ,/, 4X, 7HNBOTF- ,110) 

500 FORMAT (4F10.2) 

510 FORMAT (4X,5HRPM- , E12 . 5 , /, 4X, 5HTSS- , E12 . 5 ,/, 4X, 6HTINF- ,E12.5,/, 
14X, 6HHINF- ,E12.5) 

520 FORMAT (G10.9) 

530 FORMAT (4X,6HCRIT= ,E12.5) 

540 FORMAT (2G10.5) 

550 FORMAT ( 4X , 7HFANGL- , E12 . 5 ,/, 4X, 6HNFIN- ,15) 

560 FORMAT ( 1 5 ) 

570 FORMAT (4X,5HIFF= , 110) 

580 FORMAT (1615) 

590 FORMAT (6l5) 

600 FORMAT ( ///5X , 4HTIB= , I 5 , 10X , 3HTT- , I 5 , / , 5X , 3HTE- , I 5 , / 

1 , 6X, 3HTB-, 15 ) 

610 FORMAT (//10X, 17HCRASH , CRASH , CRASH ) 

620 FORMAT ( //5X, 4 ( E12 . 7 , 3X) ) 

630 FORMAT ( 2X , 7HSTATION, 2X , 4HQFIN , 17X , 7HQTROUGH , 1 5X , 6HQTOTAL ) 

640 FORMAT ( 4X , 1 5 , El 2 . 5 , 10X , El 2 . 5 , 10X , El 2 . 5 ) 

650 FORMAT ( // , 4X , 11HQFIN TOTAL- , E12 . 5 , 10X, 15HQTROUGH TOTAL- ,El2.5) 
660 FORMAT (/////, 4X, 5HCLI- , El 2 . 5 , 5X , 7HCANGL- , E12 . 5 ,/, 4X, 8HRBASEI- , 
1E12 . 5, 2X, 5HR2I- ,E12.5,/,4X, 8HTHICKI- , E12 . 5 , 2X , 6HBFIN- ,El2.5,/,4 
2X , 5HRPM- , E12 . 5 , 5X, 5HTSS- , E12 . 5 ,/ , 4X , 6HTINF- , E12 . 5 , 4X, 6HHINF- ,E 
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312.5,/,4X,6HCRIT= , El 2 . 5 , 4X , 7HFANGL* , E12 . 5 ,/, 4X, 5HZ0A- ,E12 
45HIFF* ,110) 

661 FORMAT ( 4X, 'FIN AREA/SMOOTH AREA* ' , El2 . 5 ) 

670 FORMAT ( 1H1 ,//2X, 4HBOA* , G12 . 5 , 5X, 4HZOA* , G12 . 5 , 5X, 5HNFIN* , 1 5 , 
1/, 5HBVIN* , G12 . 5 , 5X, 1 3HSURFACE AREA* , G12 . 5 ) 

680 FORMAT ( //5X , 2HNP , 6X , 1HX , 12X , 1HY , 12X , 1HT , 1 2X , 2HTC ) 

690 FORMAT (/2X, 13 , 3X, 4 ( FlO .6 , 3X) ) 

700 FORMAT ( //2X , 2HEL , 8X , 1HH , 11X , 9 HEL- LENGTH , 15X , 4HQ-EL ) 

710 FORMAT ( /2X , I 2 , 3X , El 2 . 5 , 3X , El 2 . 5 , 1 OX , El 2 . 5 ) 

720 FORMAT ( /2X, 22HCONVERGENCE CRI TERIAN= , El 5 . 8 ) 

730 FORMAT ( 1H , // , 5X , 4HHFG* , E12 . 5 , / , 5X , 11HNO . OF FINS* , I 5 , / , 5X , 
16HH-OUT* , El 2 . 5 , / , 5X , 5HTSAT* , El 2 . 5 , / , 5X , 4HRPM* , El 2 . 5 , / , 5X , 6HQ 
2E12 . 5 ,/, 5X, 6HQFIN = , E12 . 5 ,/, 5X, 11HHALF-ANGLE- , F8 . 3 ) 

734 FORMAT ( /5X , ' ROOTS : ' , 5X , ' REAL PARTS ' , 15X, ' IMAGINARY PARTS') 

735 FORMAT (15X,E12.5,15X,E12.5) 

740 FORMAT ( 1H0 , 6X, 1HJ, 4X, 14HFILM THICKNESS , 8X, 2HQB , 10X, 8HMASS-T 
1/, 4X, 3HTIB , 8X , 2HTT , 10X, 2HTE , 8X, 2HTB ) 

750 FORMAT ( lH , 4X, 14 , 4X, F12 . 10 , 4X, FlO . 4 , 6X, F9 . 5 , 6X, F5 . 1 ,/, 
16X,F5.1,6X,F5.1,6X,F5.1) 

760 FORMAT ( 1H0 , 6X , 1HJ , 6X , 6HK-WALL , 4X , 6HK-FI LM , 3X , 7HDENS ITY , 4X , 9 
1FILM,/, 6X, 7HEPSILON, 5X, 6HRADIUS , 5X, 4HQINC ) 

770 FORMAT ( 1H , 4X , 1 4 , 4X , F7 . 3 , 4X , F6 . 4 , 4X, F6 . 3 , 4X , F9 . 7 , 
1/,4X,F9.7,4X,F7.5,5X,F5.1,1X,F7.3) 

END 



THIS SUBROUTINE DEFINES THE POSITIONS OF SYSTEM COORDINATE P 

SUBROUTINE COORD 

COMMON/AD S/D F( 21 ) ,G(10) , IDG(100) , I GRAD , INFO , I OPT, IONED , IPRIN 
: I STRAT , IWK (2000) ,IZ(30) ,OBJ,S(2) ,VLB(2) ,VUB(2) ,W(21,30),WK(5 
:NCOLA, NCON, ND V , NGT, NRA , NR TWK , NRWK 

COMMON/OLLIE/A( 200, 50) ,AMTOT(200) ,APS,B(3) , BFIN , BOA, BVIN, C ( 3 
: CANGL , CF ( 200 ) ,CK,CLI,COF( 5) ,CW(200) , DEL (200) ,DMDOT(200) ,EA( 3 
: EPS (200) , EZERO, F( 200,1 ) , FANGL , FNOB J ( 100 ) ,H(200) , HINF , HZ ( 200 ) 
:QB( 200 ) , QINC ( 200 ) ,QTINC(200) , QTOT , QTOTAL (100) , R ( 200 ) , RB ( 200 ) 
: RBASEI , R2I , RHOF (200 ) ,ROOTI ( 4 ) ,RPM,ROOTR( 4 ) ,T( 200 ) , TALFA , TB ( 2 
:TCC( 200 ) , TE ( 200 ) , THICK , THICKI , TIB ( 200 ) , TINF , TS ( 200 ) ,TSAT,TSS 
:TT( 200 ) ,TZ ,UF( 200 ) ,X( 200 ) ,XCOF( 5) ,XPLOT( 200),Y(200),Z(200),Z 
: DOBF, DOTH, ICOR( 200,3 ) , IFF, JINT , JLC , JTC , KFF ( 50) ,KFIN( 50) ,KT,N 
: NEL , NFIN , NSNP 



C DELH IS THE STANDARD DIVISION OF FIN HEIGHT 

DELH*S( 1 )/( 12*DOBF ) 

X(l)-0.0 

Y(1)-THICK+(S(1)/12) 

N-l 

DO 20 1*1, IFF 
ICA*KFIN( I ) 

ICB*KFF( I ) 

CBA* FLOAT ( ICB-ICA) 

AN=0 .0 

DO 10 I I=ICA, ICB 

X( II )=X( 1 ) +N*AN*DELH*TALFA/CBA 
Y ( I I ) =Y( 1 ) -N*DELH 
10 AN=AN+1 . 0 
20 N=N+ 1 
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AN-0.0 

ICD-ICB-ICA+1 

DO 50 J-JTC, JLC, JINT 

X( J)-X(l) 

Y( J ) » ( 1 . O-AN/DOTH ) *TH I CK 
DO 30 J J = 1 , I CD 

X( J+JJ)-X( J )+JJ*EZERO/( 2* ( CBA+1 . 0 ) ) 
30 Y(J+JJ)=Y(J) 

JJ-ICD 

DO 40 K-1,KT 

X( J+JJ+K)=X( J+J J ) +K*APS/ ( 2 . 0*KT) 

40 Y( J+JJ+K)«Y( J) 

50 AN-AN+1.0 
RETURN 
END 



THIS SUBROUTINUE IS USED TO FORMULATE THE FEM EQUATIONS 

SUBROUTINE FORMAF 

COMMON/ADS/D F ( 21 ) ,G(10) ,IDG(100) , I GRAD , INFO , I OPT , IONED , I PRINT , 

: ISTRAT, IWK( 2000) , IZ( 30) ,OBJ,S( 2) ,VLB( 2 ) ,VUB( 2 ),W(21,30),WK(5000), 
: NCOLA , NCON , NDV , NGT , NRA , NRIWK , NRWK 

COMMON/OLLIE/A( 200,50) ,AMTOT(200) ,APS,B(3) , BFIN , BOA , BVIN , C ( 3 ) , 

: CANGL , CF ( 200 ) , CK , CLI , COF ( 5 ) ,CW(200) , DEL (200) ,DMDOT(200) ,EA(3,3) , 

: EPS (200) ,EZERO,F( 200,1) , FANGL, FNOBJ ( 100 ) ,H(200) , HINF , HZ ( 200 ) 

:QB ( 200 ) ,QINC( 200) ,QTINC( 200) , QTOT , QTOTAL (100) , R( 200 ) , RB ( 200 ) , 
:RBASEI,R2l,RHOF( 200) , ROOTI ( 4 ) , RPM , ROOTR( 4 ) ,T(200) , TALFA, TB ( 2 00 ) , 
:TCC( 200 ) , TE ( 200 ) , THICK, THICKI , TIB ( 200 ) , TINF , TS ( 200 ) ,TSAT,TSS,. 

:TT( 200) ,TZ,UF( 200) ,X( 200) ,XCOF( 5) ,XPLOT( 200) , Y( 200 ) , Z ( 200 ) ,ZOA, 
:DOBF , DOTH, ICOR( 200,3) , IFF , JINT , JLC , JTC , KFF ( 50 ) , KFIN ( 50) , KT , NBAN , 

: NEL , NFIN , NSNP 



DO 20 N«1,NSNP 
F(N, 1 )-0 . 0 
DO 10 MA-1 , NBAN 
10 A(N,MA)-0.0 
20 CONTINUE 

DO 70 IEL-1,NEL 
IA=-ICOR( IEL, 1 ) 

IB-ICOR( IEL, 2 ) 

IC=*ICOR( IEL, 3 ) 

B( 1 )-Y( IB ) -Y ( IC) 

B( 2)-Y( IC)-Y( IA) 

B( 3 )=-Y( IA ) - Y( IB) 

C( 1 )-X( IC)-X( IB) 

C( 2 )-X( IA)-X( IC) 

C( 3 )-X( IB ) -X( IA) 

LENGTH BETWEEN ELEMENT NODES 1 AND 2 
EL-SQRT(C( 3)**2+B( 3)**2) 

AREA OF TRIANGULAR ELEMENT 
AS=ABS ( (B(l)*C(2)-B(2)*C(l))/2.0) 
HC=H ( I EL ) /CK 
DO 60 J=1 , 3 
JJ=ICOR( IEL, J ) 

DO 50 K= 1 , 3 
KK=ICOR( IEL, K) 
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C FORMING THE A MATRIX 

EA( J,K)=(B( J)*B(K)+C( J)*C(K) )/(4*AS) 

IF (HC.EQ.0.0) GO TO 40 

HEL=HC*EL/6 . 0 

IF (J.EQ.3) GO TO 40 

IF (K.EQ.3) GO TO 40 

IF (J.EQ.K) GO TO 30 

EA( J , K ) = EA( J , K ) +HEL 

GO TO 40 

30 EA( J , K ) =EA ( J , K ) + 2 *HEL 
40 IF (KK.LT.JJ) GO TO 50 
NW=KK-JJ+1 

A( JJ,NW)=A( JJ,NW)+EA( J,K) 

50 CONTINUE 
60 CONTINUE 

C FORMING THE F MATRIX 

FE=HC*TS ( IEL ) *EL/2 . 0 
F ( IA, 1 ) = F( IA , 1 ) +FE 
F(IB,1)=F(IB,1)+FE 
70 CONTINUE 
RETURN 
END 



* EQUATION SOLVER OF A SYMMETRIC MATRIX THAT HAS BEEN TRANS - 

* FORMED INTO BANDED FORM. 

* 

SUBROUTINE BANDEC ( NEQ , MAXB , NVEC ) 

* 

COMMON/ADS/DF (21) ,G(10) , IDG(100) , I GRAD , INFO , I OPT , IONED , I PR IN 
: ISTRAT, IWK( 2000 ) , IZ(30) ,OBJ,S(2) ,VLB(2) ,VUB(2) ,W( 21 , 30 ) ,WK( 5 
: NCOLA , NCON , NDV , NGT , NRA, NRIWK , NRWK 

* 

COMMON/OLLIE/A( 200, 50 ) ,AMTOT(200) ,APS,B(3) , BFIN , BOA, BVIN, C ( 3 
:CANGL,CF( 200) , CK , CLI , COF ( 5 ) ,CW( 200) , DEL ( 200 ) , DMDOT ( 200 ) , EA( 3 
: EPS ( 200 ) , EZERO , F ( 200,1 ) , FANGL , FNOBJ ( 100 ) ,H( 200 ) , HINF , HZ ( 200 ) 
:QB ( 200 ) , QINC ( 200 ) ,QTINC(200) , QTOT , QTOTAL (100) , R( 200 ) , RB ( 2 00 ) 
: RBASEI , R2I , RHOF ( 200 ) ,ROOTI(4) , RPM , ROOTR ( 4 ) ,T(200) , TALFA, TB ( 2 
: TCC ( 200 ) , TE ( 200 ) , THICK , THICKI , TIB ( 200 ) , TINF , TS ( 200 ) ,TSAT,TSS 
sTT( 200) ,TZ,UF(200) ,X(200) ,XCOF(5) ,XPLOT(200) ,Y(200),Z(200),Z 
jDOBF,DOTH, ICOR( 200,3) , I FF , JINT , JLC , JTC , KFF ( 50 ) , KFIN ( 50 ) , KT , N 
: NEL , NFIN , NSNP 



LOOP*NEQ-l 
DO 20 1=1 , LOOP 
MB-I+1 

NB=MIN0( I+MAXB— 1 , NEQ ) 

DO 20 J=*MB , NB 
L= J+2-MB 
D-A(I,L)/A(I,1) 

DO 10 MM=1 , NVEC 
10 F ( J , MM ) =F ( J , MM ) -D* F ( I , MM ) 
MM-MIN0 ( MAXB-L+1 ,NEQ-J+1 ) 
DO 20 K= 1 , MM 
NN=L+K-1 

20 A( J , K ) =A( J , K ) — D*A( I , NN ) 

DO 30 1=1, NVEC 

30 F ( NEQ , I ) = F ( NEQ , I ) /A ( NEQ , 1 ) 
DO 50 1=2, NEQ 
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J-NEQ-I+1 

K-MINO ( NEQ-J+1 , MAXB ) 

DO 50 MM-1,NVEC 
DO 40 L-2,K 
MB-J+L-1 

40 F ( J , MM ) - F ( J , MM ) -A ( J , L ) * F ( MB , MM ) 
50 F( J,MM)-F( J,MM)/A( J,l) 

RETURN 

END 



SUBROUTINE DPOLRT COMPUTES THE ROOTS OF A REAL 
POLYNOMIAL USING A NEWTON-RAPHSON ITERATIVE 
TECHNIQUE. 



SUBROUTINE DPOLRT ( M, IER) 

COMMON/ADS/DF ( 21 ) ,G(10) , IDG(IOO) , IGRAD , INFO , IOPT , IONED , I PRINT , 

: I STRAT, IWK( 2000 )i , IZ( 30) ,OBJ,S( 2 ) ,VLB( 2) ,VUB( 2 ) , W( 21 , 30 ) , WK ( 500 0 ) , 
: NCOLA , NCON , NDV , NGT , NRA , NRIWK , NRWK 

COMMON/OLLIE/A( 200,50) , AMTOT ( 200) ,APS,B( 3 ) , BFIN , BOA , BVIN , C ( 3 ) , 
:CANGL ,CF( 200) ,CK,CLI,COF( 5) , CW ( 200 ) , DEL ( 200 ) , DMDOT (200), EA (3,3) , 

: EPS (200) ,EZERO,F( 200,1) , FANGL, FNOBJ ( 100 ) ,H(200) , HINF , HZ ( 200 ) , 

:QB( 200 ) , QINC ( 200 ) ,QTINC(200) , QTOT , QTOTAL (100) ,R(200) ,RB(200) , 
:RBASEI , R2I , RHOF( 200 ) ,ROOTI(4) , RPM, ROOTR( 4 ) ,T(200) , TALFA, TB ( 20 0 ) , 
:TCC( 200 ) , TE ( 200 ) , THICK, THICKI , TIB ( 200.) , TINF, TS ( 200 ) ,TSAT,TSS, 

: TT ( 200 ) , TZ , UF (200) ,X(200) ,XCOF(5) ,XPLOT(200) ,Y(200),Z(200) ,ZOA, 

: DOBF , DOTH , ICOR( 200 , 3 ) , I FF , JINT , JLC , JTC , KFF ( 50 ) , KFIN ( 50 ) , KT , NBAN , 

: NEL , NFIN , NSNP 



IFIT-0 

N-M 

IER-0 

IF ( XCOF ( N+l ) ) 10,40,10 
10 IF (N) 20,20,60 

SET ERROR CODE TO 1 

20 IER-1 
30 RETURN 

SET ERROR CODE TO 4 

40 IER=4 
GO TO 30 

SET ERROR CODE TO 2 

50 IER-2 
GO TO 30 

60 IF ( N-36 ) 70,70,50 
70 NX=N 
NXX=*N+1 
N2=l 
KJl=N+l 
DO 80 L=-1,KJ1 
MT-KJl-L+1 
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non non non non non nno 



80 C0F(MT)-XC0F(L) 



SET INITIAL VALUES 

90 X0-. 00500101 
YO-0. 01000101 

ZERO INITIAL VALUE COUNTER 
IN-0 

100 XX-XO 

INCREMENT INITIAL VALUES AND COUNTER 

XO— 10.0*YO 
YO— 10 . 0*XX 

SET X AND Y TO CURRENT VALUE 

XX-XO 
YY-YO 
IN-IN+1 
GO TO 120 
110 IFIT-1 
XPR-XX 
YPR-YY 

EVALUATE POLYNOMIAL AND DERIVATIVES 

120 ICT-0 
130 UX-0.0 . 

UY-0.0 
V-0.0 
YT-0 .0 
XT-1.0 
U-COF ( N+l ) 

IF (U) 140,270,140 
140 DO 150 1-1 , N 
L-N-I+l 

XT2«XX*XT-YY*YT 
YT2-XX*YT+YY*XT 
U-U+COF ( L ) *XT2 
V-V+COF ( L ) * YT2 
FI-I 

UX-UX+FI *XT*COF ( L ) 

UY-UY-FI * YT*COF ( L ) 

XT-XT2 
150 YT-YT2 

SUMSQ=UX*UX+UY*UY 
IF (SUMSQ) 160,230,160 
160 DX-( V*UY-U*UX) /SUMSQ 
XX-XX+DX 

DY— ( U*UY+V*UX ) /SUMSQ 
YY-YY+DY 

IF (ABS(DY)+ABS(DX)-1.0E-05) 210,170,170 

STEP ITERATION COUNTER 

170 ICT-ICT+1 

IF ( ICT-500 ) 130,180,180 
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180 IF ( I FIT ) 210,190,210 
190 IF (IN-5) 100,200,200 

( 

SET ERROR CODE TO 3 

( 

200 IER-3 

GO TO 30 

210 DO 220 L-1,NXX 
MT-KJl-L+1 
TEMP-XCOF(MT) 

XCOF ( MT ) -COF ( L ) 

220 COF(L)-TEMP 
ITEMP-N 
N-NX 

NX- I TEMP 

IF ( IFIT) 250,110,250 
230 IF (IFIT) 240,100,240 
240 XX-XPR 
YY-YPR 
250 IFIT-0 

IF ( AB S ( YY/XX ) - 1 . 0 E- 0 4 ) 280,260,260 
260 ALPHA- XX+XX 

SUMSQ-XX*XX+YY*YY 
N-N-2 
GO TO 290 
270 XX-0.0 
NX-NX-1 
NXX-NXX-1 
280 YY-0.0 

SUMSQ-0 . 0 

ALPHA-XX 

N-N-l 

290 COF ( 2 ) -COF ( 2 ) +ALPHA*COF ( 1 ) 

DO 300 L-2 , N 

300 COF ( L+l ) -COF ( L+l ) +ALPHA*COF ( L ) -SUMS Q* COF ( L-l ) 
310 ROOTI ( N2 ) -YY 
ROOTR ( N2 ) -XX 
N2-N2+1 

IF ( SUMSQ ) 320,330,320 
320 YY— YY 

SUMSQ-0. 0 
GO TO 310 

330 IF (N) 30,30,90 
END 
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